#include <feel/options.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.hpp>
#include <feel/feeldiscr/region.hpp>
#include <feel/feelpoly/im.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelpoly/polynomialset.hpp>
#include <feel/feelvf/vf.hpp>


#include <feel/feelfilters/gmshtensorizeddomain.hpp>
#include <feel/feelfilters/exporterquick.hpp>

#include <feel/feeldiscr/operatorlagrangep1.hpp>

#include <feel/feelfilters/geotool.hpp>

//include "../Schwarz2/loadmesh.hpp"

/** use Feel namespace */
using namespace Feel;
using namespace Feel::vf;


inline
po::options_description
makeOptions()
{
  po::options_description code1doptions("Code1D options");
  code1doptions.add_options()
    ("hsize1", po::value<double>()->default_value( 0.1 ), "mesh size1")
    ("hsize2", po::value<double>()->default_value( 0.1 ), "mesh size2")
    ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh( 2.0 or 2.1)")
    ;
  return code1doptions.add( Feel::feel_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Feel::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
  AboutData about( "Test_Order2Geo" ,
		   "Test_Order2Geo" ,
		   "0.2",
		   "Code 1D for -Delta u = f",
		   Feel::AboutData::License_GPL,
		   "Copyright (c) 2010 Universite Joseph Fourier");

  about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
  return about;

}


template<uint Dim,uint N,uint OrderGeo>
class Localisation
  :
  public Application
{
  typedef Application super;
public:

  //! Polynomial order \f$P_2\f$
  //    static const uint16_type Order = 5;
  ///static const uint Order = N;

  typedef double value_type;

  typedef Backend<value_type> backend_type;
  typedef boost::shared_ptr<backend_type> backend_ptrtype;

  typedef backend_type::sparse_matrix_type sparse_matrix_type;
  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
  typedef backend_type::vector_type vector_type;
  typedef backend_type::vector_ptrtype vector_ptrtype;

  typedef Simplex<Dim,OrderGeo,Dim> convex_type;
  //typedef SimplexProduct<Dim,OrderGeo,Dim> convex_type;
  typedef Mesh<convex_type> mesh_type;
  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;

  typedef bases<Lagrange<N,Scalar,PointSetFekete> > basis_type;
  //  typedef bases<Lagrange<N,Scalar/*,PointSetFekete*/> > basis_type;
  typedef FunctionSpace<mesh_type, basis_type> space_type;
  typedef boost::shared_ptr<space_type> space_ptrtype;
  typedef typename space_type::element_type element_type;

  typedef Exporter<mesh_type> export_type;
  typedef boost::shared_ptr<export_type> export_ptrtype;

  /**
   * Constructor
   */
  Localisation( int argc, char** argv, AboutData const& ad,
		po::options_description const& od )
    :
    super( argc, argv, ad, od ),
    M_backend( backend_type::build( this->vm() ) ),
    //M_backend( backend_type::build( BACKEND_TRILINOS ) ),
    meshSize1( this->vm()["hsize1"].template as<double>() ),
    meshSize2( this->vm()["hsize2"].template as<double>() ),
    exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
  {
    //Exporter<mesh_type>::New( this->vm(), this->about().appName() );
  }

  mesh_ptrtype createMesh( std::string,std::string );

  mesh_ptrtype createMesh( double meshSize );

  mesh_ptrtype createMesh2( double meshSize );
    
  mesh_ptrtype createMeshPPP( double meshSize );

  void run();
  void run2();
  void run3();


private:

  void testEvaluate(boost::tuple<mesh_ptrtype,mesh_ptrtype> __mesh);
  void testAssemblage(boost::tuple<mesh_ptrtype,mesh_ptrtype> __mesh);

  void exportResults( element_type& u,std::string __str = "u" );

private:

  // linear algebra backend
  backend_ptrtype M_backend;

  // mesh characteristic size
  double meshSize1;
  double meshSize2;

  // exporter factory
  export_ptrtype exporter;
}; // Interpolate

template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMeshPPP( double hsize )
{
    Gmsh __gmsh;
    std::string fname;
    std::ostringstream ostr;
    std::ostringstream nameStr;


    //typename imesh<Dim, Order, RDim>::ptrtype mesh( new typename imesh<Dim, Order, RDim>::type );
    mesh_ptrtype mesh( new mesh_type );


    GmshTensorizedDomain<Dim,OrderGeo,Dim,Simplex> td;
    td.setCharacteristicLength( hsize );
    td.setX( std::make_pair( -1, 1 ) );
    td.setY( std::make_pair( -1, 1 ) );

    td.setOrder( OrderGeo );
    fname = td.generate( Simplex<Dim,OrderGeo,Dim>::name() );

    ImporterGmsh< mesh_type> import( fname );
    mesh->accept( import );

    mesh->components().set( MESH_RENUMBER | MESH_UPDATE_FACES | MESH_UPDATE_EDGES );
    mesh->updateForUse();
    return mesh;
}


template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( std::string __name,std::string __str) 
{
  mesh_ptrtype mesh( new mesh_type );

  Gmsh gmsh;
  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( __name, __str  );

  ImporterGmsh<mesh_type> import( fname );
    
  float gmsh_version = this->vm()["gmsh"].template as<float>();
  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );

  return mesh;

}


template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( double meshSize ) {
  mesh_ptrtype mesh( new mesh_type );

  std::ostringstream ostr;
    
  GeoTool::Node x1(-1,-1);
  GeoTool::Node x2(0.2,1);
  GeoTool::Rectangle R1(0.1,"LEFT",x1,x2);
  R1.setMarker(_type="line",_name="Interface",_marker2=true);
  R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R1.setMarker(_type="surface",_name="Omega1",_marker1=true);
   
  Gmsh gmsh;

  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( "domain", R1.geoStr()  );

  ImporterGmsh<mesh_type> import( fname );
    
  float gmsh_version = this->vm()["gmsh"].template as<float>();

  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );

  return mesh;

}



template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::exportResults( element_type& U, std::string __str )
{

  export_ptrtype exporterLoc( Exporter<mesh_type>::New( this->vm(), this->about().appName()+__str ) );

  Log() << "exportResults starts\n";
  exporterLoc->step(0)->setMesh( U.functionSpace()->mesh() );
  exporterLoc->step(0)->add( __str, U );
  exporterLoc->save();
  Log() << "exportResults finish\n";
} // Interpolate::export


template<uint Dim,uint N,uint OrderGeo>
void
Localisation<Dim,N,OrderGeo>::testEvaluate(boost::tuple<mesh_ptrtype,mesh_ptrtype> __mesh )
{

  mesh_ptrtype mesh1 = boost::get<0>(__mesh);
  mesh_ptrtype mesh2 = boost::get<1>(__mesh);

  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh2 = space_type::New( mesh2 );

  element_type u1( Xh1, "u1" );
  element_type u2( Xh2, "u2" );


  AUTO (e ,exp(Px()*Py())*sin(2*M_PI*Px()));
  AUTO (f , sin(2*M_PI*Px())*sin(2*M_PI*Py()));
  //AUTO (g , sin(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
  //AUTO (g , Px()*(Px()-1)*Py()*(Py()-1) );

  //AUTO (h , Px()*(Py()+1) );
  //AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );
  AUTO (h , cst(1.)*Px());

  AUTO (g , 2.5*Px()+0.5);

  u1 = project( Xh1, elements(mesh1), h );
  u2 = project( Xh2, elements(mesh2), h );

  //-----------------------------------------------------------------------------------//

  double __errId = 
    std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")), idv(u1)*idv(u1) ).evaluate()(0,0))
    -std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")), idv(u1)*idv(u1) ).evaluate()(0,0));

  double __errGrad = 
    std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")), gradv(u1)*trans(gradv(u1)) ).evaluate()(0,0))
    -std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")), gradv(u1)*trans(gradv(u1)) ).evaluate()(0,0));

  double __errId2 = 
    std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")), idv(u2)*idv(u2) ).evaluate()(0,0))
    -std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")), idv(u2)*idv(u2) ).evaluate()(0,0));

  double __errGrad2 = 
    std::sqrt(integrate( markedfaces(mesh1, mesh1->markerName("Interface")), gradv(u2)*trans(gradv(u2)) ).evaluate()(0,0))
    -std::sqrt(integrate( markedfaces(mesh2, mesh2->markerName("Interface")), gradv(u2)*trans(gradv(u2)) ).evaluate()(0,0));


  std::cout << "\n -----u1 id boundary Interface :" 
	    << __errId << "\n";

  std::cout << "\n -----u1 grad boundary Interface :" 
	    << __errGrad << "\n";

  std::cout << "\n -----u2 id boundary Interface :" 
	    << __errId2 << "\n";

  std::cout << "\n -----u2 grad boundary Interface :" 
	    << __errGrad2 << "\n";

}

template<uint Dim,uint OrderChamp,uint OrderGeo>
void
Localisation<Dim,OrderChamp,OrderGeo>::run()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

  using namespace Feel::vf;  
  //-----------------------------------------------------------------------------------//

  int MyOrder=OrderChamp;

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % MyOrder
			  % this->vm()["hsize1"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//

  //Geometry
  GeoTool::Node x1(-1,-1);
  GeoTool::Node x2(0,1);
  GeoTool::Rectangle R1( meshSize1,"LEFT",x1,x2);
  R1.setMarker(_type="line",_name="Interface",_marker2=true);
  R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R1.setMarker(_type="surface",_name="Omega1",_marker1=true);
 
  GeoTool::Node x3(0,-1);
  GeoTool::Node x4(1,1);
  GeoTool::Rectangle R2(meshSize2,"RIGHT",x3,x4);
  R2.setMarker(_type="line",_name="Interface",_marker4=true);
  R2.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true);
  R2.setMarker(_type="surface",_name="Omega2",_marker1=true);

  //-----------------------------------------------------------------------------------//

  //Mesh
  mesh_ptrtype mesh1 = createMesh( "mesh1",R1.geoStr() );
  mesh_ptrtype mesh2 = createMesh( "mesh2",R2.geoStr() );


  std::list<boost::tuple<mesh_ptrtype,mesh_ptrtype> > __listMesh;
  __listMesh.push_back(boost::make_tuple(mesh1,mesh2));

  typename std::list<boost::tuple<mesh_ptrtype,mesh_ptrtype> >::iterator itMesh = __listMesh.begin();
  typename std::list<boost::tuple<mesh_ptrtype,mesh_ptrtype> >::iterator itMesh_end = __listMesh.end();

  for ( ; itMesh!=itMesh_end ; ++itMesh)
    {
      testEvaluate( *itMesh);
    }


  std::cout << "\n---------------------------------------"
	    << "\n---------------------------------------"
	    << "\n---------------------------------------\n";


  mesh_ptrtype mesh1bis = createMesh( "mesh1bis",R1.geoStr() );
  space_ptrtype Xh1 = space_type::New( mesh1 );
  space_ptrtype Xh1bis = space_type::New( mesh1bis );

  element_type u1( Xh1, "u1" );
  element_type u1bis( Xh1bis, "u1bis" );

  AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );
  //AUTO (h ,exp(Px()*Py())*sin(2*M_PI*Px()));

  //AUTO (h , cst(1.) );


  u1 = project( Xh1, elements(mesh1), h );
  u1bis = project( Xh1bis, elements(mesh1bis), h );

  /*
  double __errId = 
    std::sqrt(integrate( elements(mesh1), idv(u1)*idv(u1) ).evaluate()(0,0))
    -std::sqrt(integrate( elements(mesh1bis), idv(u1)*idv(u1) ).evaluate()(0,0));

  double __errGrad = 
    std::sqrt(integrate( elements(mesh1), gradv(u1)*trans(gradv(u1)) ).evaluate()(0,0))
    -std::sqrt(integrate( elements(mesh1bis), gradv(u1)*trans(gradv(u1)) ).evaluate()(0,0));

  std::cout << "\n -----u1 id boundary Interface :" 
	    << __errId << "\n";

  std::cout << "\n -----u1 grad boundary Interface :" 
	    << __errGrad << "\n";
  */

  /*
  const uint visOrder=3;
  typedef bases<Lagrange< visOrder, Scalar, PointSetFekete> > vis_basis_type;
  typedef FunctionSpace< mesh_type, vis_basis_type,Continuous, double> vis_functionspace_type;
  typedef boost::shared_ptr<vis_functionspace_type> vis_functionspace_ptrtype;
  typedef typename vis_functionspace_type::element_type vis_element_type;
  
  typedef OperatorLagrangeP1<vis_functionspace_type> operator_lagp1_type;
  typedef typename boost::shared_ptr<operator_lagp1_type> operator_lagp1_ptrtype;
  
  typedef typename operator_lagp1_type::dual_image_space_ptrtype vis_space_ptrtype;

  
  vis_functionspace_ptrtype Wh;
  operator_lagp1_ptrtype op_lagp1;
  vis_space_ptrtype vis_space;

  Wh = vis_functionspace_type::New( Xh1->mesh() );
  op_lagp1 = operator_lagp1_ptrtype(new operator_lagp1_type( Wh, M_backend ));
  */
  //vis_space = op_lagp1->dualImageSpace();
  
  //vis_element_type = project(vis_space

  //vis_space->mesh();



  /*  mesh_ptrtype meshhh = createMeshPPP(0.1);
  boost::shared_ptr<space_type> Xh( new space_type(meshhh) );
  typename space_type::element_type u( Xh, "ooioi" );
  u = project( Xh, elements(meshhh), Px() );
  */

  typedef Backend<value_type> backend_type;
  boost::shared_ptr<backend_type> backend( backend_type::build( BACKEND_GMM ) );

  OperatorLagrangeP1<space_type> I( Xh1, backend );
  typedef typename OperatorLagrangeP1<space_type>::dual_image_space_type::mesh_type image_mesh_type;
  typename OperatorLagrangeP1<space_type>::dual_image_space_ptrtype Yh( I.dualImageSpace() );
  typename OperatorLagrangeP1<space_type>::dual_image_space_type::element_type yy( Yh, "yy_elem" );
  yy = project( Yh, elements(Yh->mesh()), Px() );

  /*
  export_ptrtype exporterLoc( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) );
  exporterLoc->step(0)->setMesh( vis_space->mesh() );
  exporterLoc->step(0)->add( "oollaa", u1 );
  exporterLoc->save();
  */
  //ExporterQuick<image_mesh_type> exp( "olol", "ensight" );
  //exp.save( 0,  yy );

  //  exportResults(yy,"yu1");
  //exportResults(u1,"u1");

}


template<uint Dim,uint OrderChamp,uint OrderGeo>
void
Localisation<Dim,OrderChamp,OrderGeo>::run2()
{


  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_type::name()
			  % OrderChamp
			  % this->vm()["hsize1"].template as<double>()
			  );


  //Geometry
  GeoTool::Node x1(0,0);
  GeoTool::Node x2(1,1);
  GeoTool::Rectangle R1( meshSize1,"LEFT",x1,x2);
  R1.setMarker(_type="line",_name="Interface",_marker2=true);
  R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R1.setMarker(_type="surface",_name="Omega1",_marker1=true);

  GeoTool::Node x_centre(0.5,0.5);
  GeoTool::Node x_bord(0.8,0.5);
  GeoTool::Circle C(meshSize2,"uncercle",x_centre,x_bord);
  C.setMarker(_type="line",_name="trou",_marker1=true);

  //Mesh
  mesh_ptrtype mesh = createMesh( "mesh1",(R1-C).geoStr() );
  //mesh_ptrtype mesh = createMeshPPP(meshSize1);

  //geomap
  typedef typename mesh_type::gm_type gm_type;
  typedef typename boost::shared_ptr<gm_type> gm_ptrtype;
  gm_ptrtype __gm = mesh->gm();
 
  //geomap context
  typedef typename mesh_type::element_type geoelement_type;
  typedef typename gm_type::template Context<vm::POINT, geoelement_type> gmc_type;
  typedef boost::shared_ptr<gmc_type> gmc_ptrtype;

  //pts de Gauss
  typename _Q<4>::template apply<mesh_type::nDim,typename mesh_type::value_type, Simplex >::type im;

  // iterators sur les elements du maillage
  typename mesh_type::element_iterator el_it;
  typename mesh_type::element_iterator el_en;
  boost::tie( boost::tuples::ignore, el_it, el_en ) = Feel::elements( *mesh );

  typedef typename mesh_type::gm_type::precompute_type geopc_type;
  typedef typename mesh_type::gm_type::precompute_ptrtype geopc_ptrtype;
  geopc_ptrtype __geopc( new geopc_type( __gm, im.points() ) );

  for ( ;el_it!=el_en;++el_it) 
    {
      //init and apply the geometric transformation phi
      gmc_ptrtype __c( new gmc_type( __gm,
				     *el_it,//mesh->element( *el_it ),
				     __geopc ) );

      //init inverse geometric transformation phi^{-1}
      typename mesh_type::Inverse::gic_type gic( __gm, *el_it );

      for (uint i=0;i< __c->xReal().size2();++i)
	{

	  // get phi for one point
	  typename mesh_type::node_type n = ublas::column( __c->xReal(), i );

	  //compute phi^{-1}
	  gic.setXReal(n);

	  //verify that : phi°phi^{-1}=Id
	  for (uint j=0;j<n.size();++j)
	    {
	      double err = std::abs(ublas::column(im.points(),i)(j) - gic.xRef()(j));
	      if ( err> 1e-9) std::cout<< "\nProb : "<< err;
	    }

	}

    }
  std::cout<<"\nfinish test geomap inverse!\n";


  space_ptrtype Xh = space_type::New( mesh );

  element_type u( Xh, "u" );

  //AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );
  AUTO (h ,exp(Px()*Py())*sin(2*M_PI*Px()));

  u = project( Xh, elements(mesh), h );

  exportResults(u,"u");


}

template<uint Dim,uint OrderChamp,uint OrderGeo>
void
Localisation<Dim,OrderChamp,OrderGeo>::run3()
{
  //Geometry
  GeoTool::Node x1(0,0);
  GeoTool::Node x2(1,1);
  GeoTool::Rectangle R1( meshSize1,"LEFT",x1,x2);
  R1.setMarker(_type="line",_name="Interface",_marker2=true);
  R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
  R1.setMarker(_type="surface",_name="Omega1",_marker1=true);

  GeoTool::Node x_centre(0.5,0.5);
  GeoTool::Node x_bord(0.8,0.5);
  GeoTool::Circle C(meshSize2,"uncercle",x_centre,x_bord);
  C.setMarker(_type="line",_name="trou",_marker1=true);

  std::ostringstream ostr;
  ostr << "Mesh.MshFileVersion = 2;\n"
       << "h=" << meshSize1 << ";\n"
       << "Point(1) = {0,0,0, h};\n"
       << "Point(2) = {1,0,0, h};\n"
       << "Point(3) = {1,1,0, h};\n"
       << "Point(4) = {0,1,0, h};\n"
       << "Line(1) = {4, 3};\n"
       << "Line(2) = {3, 2};\n"
       << "Line(3) = {2, 1};\n"
       << "Line(4) = {1, 4};\n"
       << "Line Loop(9) = {8, 5, 6, 7};\n"
       << "Line Loop(10) = {1, 2, 3, 4};\n"
       << "Plane Surface(11) = {10};\n"
       << "Extrude {0, 0, 1} {\n"
       << "Surface{11};\n}"
       << "Surface Loop(34) = {24, 11, 20, 33, 28, 32};\n"
       << "Volume(35) = {34};\n";


  //Mesh
  mesh_ptrtype mesh = createMesh( "mesh1",(R1).geoStr() );
  //mesh_ptrtype mesh = createMesh( "mesh3D",ostr.str() );

  mesh->updateForUse();
  mesh->tool_localization()->updateForUse();
  mesh->tool_localization()->kdtree().nbNearNeighbor(3);
  mesh->tool_localization()->kdtree().writeLatexData();
}


/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{
  //dim, N, OrderGeo
  Localisation<2,6,1> moncode( argc, argv, makeAbout(), makeOptions() );

  //moncode.run();
  //moncode.run2();
  moncode.run3();
}



